Tracking global invasion pathways of the spongy moth (Lepidoptera: Erebidae) to the United States using stable isotopes as endogenous biomarkers

Abstract The spread of invasive insect species causes enormous ecological damage and economic losses worldwide. A reliable method that tracks back an invaded insect's origin would be of great use to entomologists, phytopathologists, and pest managers. The spongy moth (Lymantria dispar, Linnaeus 1758) is a persistent invasive pest in the Northeastern United States and periodically causes major defoliations in temperate forests. We analyzed field‐captured (Europe, Asia, United States) and laboratory‐reared L. dispar specimens for their natal isotopic hydrogen and nitrogen signatures imprinted in their biological tissues (δ2H and δ15N) and compared these values to the long‐term mean δ2H of regional precipitation (Global Network of Isotopes in Precipitation) and δ15N of regional plants at the capture site. We established the percentage of hydrogen–deuterium exchange for L. dispar tissue (Pex = 8.2%) using the comparative equilibration method and two‐source mixing models, which allowed the extraction of the moth's natal δ2H value. We confirmed that the natal δ2H and δ15N values of our specimens are related to the environmental signatures at their geographic origins. With our regression models, we were able to isolate potentially invasive individuals and give estimations of their geographic origin. To enable the application of these methods on eggs, we established an egg‐to‐adult fraction factor for L. dispar (Δegg‐adult = 16.3 ± 4.3‰). Our models suggested that around 25% of the field‐captured spongy moths worldwide were not native in the investigated capture sites. East Asia was the most frequently identified location of probable origin. Furthermore, our data suggested that eggs found on cargo ships in the United States harbors in Alaska, California, and Louisiana most probably originated from Asian L. dispar in East Russia. These findings show that stable isotope biomarkers give a unique insight into invasive insect species pathways, and thus, can be an effective tool to monitor the spread of insect pest epidemics.


| INTRODUC TI ON
The global economic and ecological burden posed by invasive insect species is high, with costs rising to $70 billion dollars per year (Bradshaw et al., 2016). The ability to identify the natal origin of a pest insect captured on non-native territory would significantly support the implementation of pest management strategies. Therefore, the determination of invasive insects' natal origins has been a topic of investigation for many years, but has largely focused on genetic analysis (Picq et al., 2018;Wu et al., 2015Wu et al., , 2020. While genetic analyses can yield insight into a population's geographic origin, it does not provide the nuanced data often required for effective management of a newly arrived invasive species, such as whether the individual was born on site (domestic/local) or transported in (exotic). On the other hand, stable isotope analysis (SIA), particularly of hydrogen, nitrogen, sulfur, and oxygen, is a useful tool for tracking an individual's recent whereabouts (Bowen & West, 2008;Hobson, 1999).
The European spongy moth Lymantria dispar dispar (Lepidoptera: Erebidae, Linnaeus 1758) (subsequently referred to European L. dispar and formerly known as the gypsy moth) is native to all of temperate Eurasia and northern Africa, and was brought to the United States between 1868 and 1869 near Boston, Massachusetts, from where it escaped and spread shortly after (Fernald & Forbush, 1896).
In the United States, there are 400 different tree species European L. dispar larvae feed on (US Department of Agriculture et al., 1981).
Additionally, the moth faces no major natural enemies, competitors, or diseases in North America that would control populations like in their native habitat (Liebhold et al., 1995). Consequently, the moth has been able to gradually expand its distribution. The area infested by European L. dispar is confined to the northeastern United States and the eastern maritime provinces of Canada, with an advancing front slowly moving in a southwesterly direction (Epanchin-Niell & Liebhold, 2015). The population periodically builds to outbreak levels that can result in serious economic, environmental, and public nuisance problems (Liebhold et al., 2000;Thorpe et al., 2006). Since 1924, tens of millions of hectares have been defoliated in the United States forests by the European L. dispar (Sharov et al., 2002), leaving regional economies with costly pest mitigation and prevention measures (Bigsby et al., 2014;Epanchin-Niell & Liebhold, 2015;Jardine & Sanchirico, 2018).
The natural expansion of European L. dispar in North America advances relatively slowly with 13 miles/year, which primarily is a consequence of the female's flight incapability (Reineke & Zebitz, 1998). This is an advantage for pest management, as egg masses can be found very close to the pupation sites (Fernald & Forbush, 1896).
Human activities near the pupation sites, however, have significantly accelerated the spread throughout the country, either through moths placing eggs on vehicles, cargo, and gear or through displacement of infested natural resources like wood (Bigsby et al., 2011;Continental Dialogue, 2019). In addition to that, the intensifying container shipping from East Asia, including eastern China, Japan, the Republic of Korea, and eastern Russia, has opened a substantial entry pathway to North America for the Asian spongy moth Lymantria dispar asiatica (which include L. dispar asiatica [Vnukovskij] and L. dispar japonica [Motschulsky]) (subsequently referred to as Asian L. dispar) (Gray, 2017;Paini et al., 2018). One major difference between the Asian and European L. dispar is that females of the former subspecies are capable of strong and directed flight. The brightly lit shipping ports have been shown to attract Asian L. dispar females that can lay eggs on vessel superstructure and shipping containers (Gray, 2017;Wallner et al., 1995). If the Asian L. dispar were to establish in North America, they can hybridize with European L. dispar and produce fertile offspring with flight or gliding capable females, which might accelerate population expansion of the pest, as some researchers believe (Gray, 2017;Keena et al., 2007;Robinet & Liebhold, 2009). The Asian L. dispar is also adapted to colder climates and higher altitudes and therefore has an even broader host range than European L. dispar, with 500 host species, including several coniferous trees (Trotter et al., 2020). The Asian L. dispar is not yet established in the United States, but adult Asian L. dispar males of unknown natal origin were captured in several US states in the past years, which gave reasons for concern for natural biospheres and called for action to intensi- There are L. dispar surveillance and eradication programs, both inside and outside of the currently infested regions of North America (Continental Dialogue, 2019;Epanchin-Niell & Liebhold, 2015;Liebhold et al., 2021;Sharov et al., 2002;Sharov & Liebhold, 1998 (Bigsby et al., 2011;Paini et al., 2018). To effectively slow the spread, limit the establishment of L. dispar across North America and to locally eradicate new outbreaks, it will be important to determine whether newly detected moths or eggs are from a locally established population or are recent arrivals from imported cargo (Gray, 2017). In some instances, genetic analysis of a specimen may be capable of determining a linage's geographic descent (Wu, 2016;Wu et al., 2015Wu et al., , 2020. However, stable isotope analysis (SIA) can help control and eradication programs by adding important information on the geographic origin and introduction pathways of an individual, the insect's diet, and alternative host species, which can help clarify whether a detection is newly introduced or has been locally established (Hungate et al., 2016; International Atomic Energy Agency, 2009).

Biogeochemistry
The stable isotope signatures of an individual reflect the natural, geographically dependent, variations of isotopic signals in precipitation, soil, and vegetation of the individual's origin (Cameron, 2005).
If the location of the individual changes rapidly (in relation to the endogenous turnover rate of its tissue), it will retain its origin's isotopic profile while already living in the new destination (International Atomic Energy Agency, 2009). Since chitin structures like insect exoskeletons have a slow turnover rate, isotopic signals of the chitin tissue, such as the deuterium-hydrogen ratio (δ 2 H) and the 15 N/ 14 N ratio (δ 15 N), are suitable biomarkers to track an insect's recent migration (Hungate et al., 2016). Tables 1 and 2 show the mean annual δ 2 H values in precipitation and mean annual δ 15 N values in foliage depending on the geographic regions (Bowen & West, 2008;Terzer et al., 2013).
The L. dispar is exceptionally suitable for SIA to track the natal origins, as it only feeds in the larval stage (Drooz, 1985). The 40 days feeding period of larvae starts between May and July after hatching. In the short postmetamorphosis adulthood (6-10 days), moths do not possess a digestive system. They solely mate, lay the overwintering egg masses from which larvae hatch the following spring, and die (Drooz, 1985;US Department of Agriculture et al., 1981).
A L. dispar's isotopic composition therefore is a direct derivative of the vegetation's isotopic composition it fed on as a larva (Hungate et al., 2016), and it does not significantly change during the adult stage. As the eggs are formed inside the females' adult body, they reflect the mother's isotopic values (International Atomic Energy Agency, 2009;Montgomery, 1982) until they hatch and start feeding the next year.
A small proportion of an insect's hydrogen isotope value, however, is only loosely bound in the tissue and can be exchanged with hydrogen in the water of the ambient air (Hungate et al., 2016;Qi & Coplen, 2011). Thus, the isotopic value of tissue changes slightly (also postmortem) within a certain equilibration period, depending on the surrounding air conditions. This phenomenon is called hydrogen-deuterium exchange, and its magnitude is determined by the species-specific percentage of exchangeable hydrogen (Hungate et al., 2016;International Atomic Energy Agency, 2009;Wassenaar & Hobson, 2003 Our study aims to help develop stable isotope tools for pest pathway analysis and control programs in general, as the principles used in this study can be applied to all chitin tissues of insects (Bowen et al., 2005;Bowen & West, 2008;Hobson et al., 2004;Hungate et al., 2016;International Atomic Energy Agency, 2009;Mekki et al., 2016).

| MATERIAL S AND ME THODS
To assess a reproducible method to determine the natal provenance of possibly exotic Asian L. dispar or European L. dispar samples (as done before for other insect pests by Hungate et al., 2016), we measured field-captured and laboratory-reared moths and eggs for their deuterium and nitrogen isotopic signatures.  Note: For the control groups ("ctrl", Mongolian Asian Lymantria dispar "LDAM" and New Jersey strain of European L. dispar "NJSS", last three lines), the value for "local precipitation δ 2 H" is the local tap water value that was used for the artificial diet, written inside parentheses.

| Sampling and measurements
We analyzed samples together with reference materials (RMs) for their isotopic hydrogen (δ 2 H) and nitrogen (δ 15 N) ratio at the We applied the same methods for the control (Ctrl) group studies. Quadruplet samples of Ctrl-European L. dispar and Ctrl-Asian L. dispar were prepared and measured together for δ 2 H and δ 15 N.
A second set of Ctrl-Asian together with Ctrl-Asian eggs were prepared and measured afterward, by another person, to define the egg-to-adult offset.
We used the reference materials' expected versus measured linear regression to obtain a calibration equation. We corrected the insect values for the measurement inaccuracy using the equation where y is the measured sample value, k is the slope, and d the intercept of the calibration curve.

| Egg-to-adult conversion
Eggs generally showed a high negative offset from adult moths (up to 30‰). We established an egg-to-adult fraction factor (16.3 ± 4.3‰) that allows to approximate the mother's isotopic value. We added the offset to eggs' measured values after the correction of the measurement inaccuracy (see above) and analyzed them like adult moth samples from thereon.

| Hydrogen-deuterium comparative equilibration experiment
In order to obtain the fixed (nonexchangeable) value δ2H natal , we determined the percentage of exchangeable hydrogen (P ex = 1 − P nex ) in comparative equilibration experiments (Hungate et al., 2016;Wassenaar & Hobson, 2003) and removed the exchangeable part (δ2H Amb ) from the moth value. Two groups of five 0.2 ± 0.025 mg samples from Wisconsin moths' tarsus were separately equilibrated with two different microatmospheres for 1 week (δ2H Group1 = + 4.4 ‰ and δ2H Group2 = − 157 ‰) until the equilibrium δ2H sample = P ex * δ2H Amb + P nex * δ2H natal (Equation 3) was stable. Then the samples were dried in a vacuum oven at 60°C for 4 days to remove moisture and transferred to measurement immediately thereafter (Qi & Coplen, 2011). Analytes were not exposed to ambient air in the laboratory for longer than 2 h before measurements. With the two distinct equilibrated groups, we cal-

| Expected insect values
We approximated the expected chitin

| Outliers
Values with an offset greater than 20‰ (δ 2 H) or 2‰ (δ 15 N) between the individual natal moth value and the expected insect value (see above) were considered outliers. In addition, we considered values δ 15 N > 5‰ as outliers, according to the natural boundary in Bowen and West (2008). Outliers were automatically submitted to origin backtracking.

| Regression analysis
Based on trap location details and satellite images, we estimated the geographical coordinates of the moths' capture sites. These data were not recorded upon capture, but needed for the capture site's isotopic signal estimates. We retrieved the monthly averages of δ 2 H in precipitation at the moth's capture sites during the assumed feeding period from the Online Isotopes in Precipitation Calculator (OIPC) (Bowen, 2018;Bowen et al., 2005;IAEA and WMO, 2015;Welker, 2000). We corrected this expected δ 2 H insect values for climate change if moths were caught before 2000 (Liu et al., 2018). An evaluation for land-use change in the region would have been useful but was out of scope for this study.
Subsequently, we performed a correlation analysis between all individuals that had passed the outlier tests or were confirmed domestic (specimens from Minnesota, US) and the expected environmental signature at their collection site. Isolated values with a disproportionately big negative impact on the correlation were removed manually (seven samples for δ 2 H and four samples for δ 15 N).
All moth values were regressed into their corresponding δ 2 Hprecipitation and δ 15 N plant values by the obtained regression model. The goodness of fit was evaluated for the residuals between the regressed and the environmental values (Δ modeled-precipitation δ 2 H and Δ model-plant δ 15 N) in a two-dimensional plot.

| Origin analyses
To detect "likely exotic" colonies (as opposed to "likely domestic"), we performed a heteroscedastic two-tailed Student's t test at a 0.05 significance level (0.01 level for δ 2 H of field-collected adult moths due to high sensitivity) between the moths' isotopic values and the expected insect values. Outliers and values that failed the t test were matched with their corresponding δ 2 H or δ 15 N geographical zone (see Tables 1 and 2).
Additionally, USDA-intern genetic analysis of some samples (methods not shown here) (Bogdanowicz et al., 1993;Garner & Slavicek, 1996;Wu et al., 2020) gave external validation of some of the model-given predictions of origins.

| Natural isotopic variation and egg-to-adult conversion
An Asian L. dispar and a European L. dispar strain were reared at the USDA OTIS Laboratory's insectary and measured for their δ 2 H and δ 15 N values (Table 4). The average isotopic variation for moths raised under identical conditions (average distance between minimum and maximum) was 11.9‰ for δ 2 H and 1.3‰ for δ 15 N. The average eggto-adult offset for spongy moths was Δ egg δ 2 H = 16.30 ± 4.3‰ (not assessed for δ 15 N).
The offset between the Ctrl-European and Ctrl-Asian-1 medians was 15.9‰ for δ 2 H and 1.88‰ for δ 15 N, thus almost as low as the average natural variation for L. dispar. The maximum distance between Ctrl-European and a Ctrl-Asian was as high as 24.6‰ for δ 2 H and 3.5‰ for δ 15 N.

| Hydrogen-deuterium exchange correction
The percentage of exchangeable hydrogen (P ex ) for spongy moths was 8.2% (Table 5). All obtained P ex values were in between the plausibility thresholds of 6% < P ex < 14% (Qi & Coplen, 2011). The maximum observable isotopic shift of a L. dispar due to P ex correction was Δδ 2 H = ±1.67‰. The reference material Casein 139443 was the most comparable RM to spongy moth tissue in terms of P ex .
Additionally, all moths tested from Brown County, Wisconsin (US, n = 4) were far outside the expected range for both isotopes.
In Carlton County (Minnesota, US), London (U.K.), and Koshunai (Japan), the moths' δ 2 H separated into two pairs that could be from neighboring populations (not presented here). Only one pair from Minnesota was removed as outliers.

| Locally characteristic moth values
We obtained the natal δ 2 H and δ 15 N values by removing the percentage of exchangeable hydrogen. Generally, East Asian moths showed more δ 2 H-negative, but more δ 15 N-positive values than European moths ( Figure 1). However, there was a variegated overlap zone between European and Asian moths. The eggs found in US harbors resembled European values in their δ 15 N, but were much more δ 2 H negative, as already reported through the egg-to-adult offset.
Most moth colonies were satisfactorily close to their expected values, with 85% of δ 2 H values being categorized as "likely domestic" ( Figure 2). The δ 15 N values overlapped with the expected ones only in 52% of locations, which was expected due to the uncertainties discussed in Methods. The control groups did not overlap with the expected range, which was acceptable since the moths were not fed local vegetation.
With the applied egg-to-adult offset for δ 2 H, the eggs from the control group overlapped identically with the expected value ( Figure 2). Accordingly, this suggests that the eggs found in Juneau (Alaska) and Portland (Oregon) may have originated from a local reproducing population or from a region characterized by a similar climate and topography. For the eggs found in New Orleans (Louisiana) and Long Beach (California), the isotope data strongly suggested a geographically different origin. The δ 15 N measurements did not agree on the eggs from New Orleans (LA).

| Worldwide moth-origin regression
The worldwide regression analyses yielded a moderate association between δ 2 H moth and δ 2 H precipitation (R 2 = .548), as well as δ 15 N moth and δ 15 N plant (R 2 = .530) (Figure 3), after some moths were manually removed (see Methods).
The natal environment's isotopic values plotted against the moth's values could be divided into a European and Asian section with a narrow intermixture zone between −30‰ and − 40‰ for δ 2 H and 0‰ and 1‰ for δ 15 N (Figure 3). Most US sites resembled the East Asian precipitation signatures, except for Wisconsin.
We calculated the offset between the regressed moth values (by means of the regressions in Figure 3) and the environmental signatures Δ modeled-precipitation δ 2 H and Δ modeled-plant δ 15 N. The combined goodness of fit was moderately good (Figure 4). Most residuals were distributed in the ranges between −30‰ < δ 2 H natal < 30‰ and −5‰ < δ 15 N < 5‰. We assumed no systematic error in the model, since the residuals showed a random distribution.

| Origin analysis
Values that were earlier categorized as outliers, manually sorted out from the regression analysis, or that were classified as "likely exotic" by the Student's t test were regressed into their environmental δ 2 H precipitation and δ 15 N plant values and matched to the corresponding geographical zones (Table 6 and Figure 5).
The hydrogen and nitrogen analyses did not agree on outliers in most cases. Some samples from the United States (Chittenden and Brown) did not match the local expected values and were considered "likely exotic" by the model, though they were collected from local populations there. Their matching geographic zones (Tables 1 and 2), however, fit to their collection sites in the end. The eggs found on a ship in Alaska seemed to fit the local expected values perfectly at first, however, the δ 2 H t test classified the regressed values as "likely exotic" later. The combination of genetic analysis (Wu et al., 2020) and isotopic results suggested northeastern Russia as a possible origin. The eggs found on ships in Long Beach (CA) and New Orleans (LA) may originate from southern Russia as well. The eggs from a ship a All P ex values were below the reject threshold of 14% (Qi & Coplen, 2011) and were used to correct the tissue's hydrogen-deuterium exchange with ambient air (Equation 4).

TA B L E 5
Percentage exchangeable hydrogen (P ex ) for insect samples and reference materials (RMs) obtained from comparative equilibration experiments F I G U R E 1 Isotopic signatures of moths (squares and triangles) and *eggs (diamonds) with standard deviations. The rectangles are continental means and the triangles regional ones. Blue stands for moths found in Asia and red for moths found in Europe. Small green diamonds are the regional means and the large green diamond is the mean of all egg samples berthing in Portland (OR) were classified as "likely domestic" by both isotopes. Since there is no known established population in Oregon, and the genetic analysis identified an Asian L. dispar descent, this could suggest that we lack isotopic data from an Asian location that has very similar features to Oregon. Alternatively, there might be an establishing Asian L. dispar colony in Oregon, even though this hypothesis is very unlikely given diligent L. dispar survey work performed in this state. Finally, the outlier found in Corsica seemed to match East Asian environmental signatures like Japan's. This sample had not been genetically analyzed.

| DISCUSS ION
In this study, we analyzed field-collected spongy moths for their natal δ 2 H and δ 15 N imprinted in their biological tissues and compared these values to the long-term mean δ 2 H of precipitation and Isotopic analysis of hydrogen and nitrogen alone cannot guarantee that a moth/egg that resembled the expected value is in fact native at the investigated site, but it may be originating from any other place in the world that has very similar climatic, topographical, and vegetation features. Here, heavier isotope ratios, such as sulfur or strontium, may be useful in determining the true origin of the specimens, yielding information on the natal distance from the sea and underlying geology, respectively, however with increasing complexity and cost (Holder et al., 2014;Schmidt et al., 2005). This highlights that isotope studies are better at answering questions such as "Is the moth from here?" rather than open questions like "Where is the moth from?" Generally, the δ 2 H values were significantly more sensitive, and more reliable, than δ 15 N for determining the origin of L. dispar (Montgomery, 1982), also because better research tools are available (Bowen, 2018;IAEA and WMO, 2015). In some of the collec-  (Bowen & West, 2008) and for δ 2 H in precipitation (bottom row) (Terzer et al., 2013). The expected value was retrieved from the Online Isotopes in Precipitation Calculator (Bowen, 2018) Alternatively, the phenomenon can be caused by a different primary water source of the vegetation than precipitation, that is, snow or underlying groundwater, which do not necessarily reflect the precipitation's signature (Finger et al., 2013). In the case of Minnesota, for example (one grouping at −90.4‰, the other at −67.3‰), the underlying water table has a signature similar to the close-by Lake Superior (δ 2 H = −68‰), which is significantly lower compared to the OIPC precipitation values (δ 2 H = −46.3‰) (Foley et al., 2014).
The difference of one group of larvae feeding on vegetation that primarily takes up precipitation water and another group feeding on vegetation that takes underground water provided by a near-by surface water body would likely be detectable in L. dispar isotopic values (Montgomery, 1982). This not only demonstrates the high sensitivity of δ 2 H values, but also shows that any deviation from the expected value has to be assessed individually with attention to the vegetation's water sources.
Multi-isotope comparison of regressed moth/egg values with local environmental values has the potential to yielding more clearcut estimates of origins, but needs more precise isotopic regression models as input on one hand, and more reference locations on the other. We concede that our model's accuracy could be improved by implementing a correction for the different heights above sea level, a parameter that should be recorded during moth collection.  (Wu et al., 2020).
that international trade and transport play a significant role in the disseminating of L. dispar to the United States, and underlines the importance of current surveillance and tracking the invasion pathways for L. dispar. Finally, we suggest that the methods presented in this article could be applied to any insect species with a chitindominant body part of low turnover rate, and similar methods for avian feather and mammal hair have been applied in previous animal migration studies (Bowen et al., 2005;Bowen & West, 2008;Foley et al., 2014;Hobson, 1999;Holder et al., 2014Holder et al., , 2020Holder et al., 2015;Hungate et al., 2016;International Atomic Energy Agency, 2009;Mekki et al., 2016).
Based on the presented results, we recommend that insect trapping and surveillance, as well as control, mitigation, and inspection regimes at embarkation ports are of utmost importance to prevent uncontrollable, costly invasion. For using SIA in pest management, those factors should be considered: (1) Collection of the exact geographic coordinates of the insect capture sites, including height above sea level and where the sample was collected from (ship, cargo, plant, or other feature in the natural environment), and collection of regional vegetation samples from possible host species on the same day to track the regional δ 15 N plant background. (2)  Rebecca Nowotny-Hood: Conceptualization (lead); data curation (equal); formal analysis (equal); funding acquisition (lead); investigation (supporting); methodology (lead); project administration (lead); resources (lead); supervision (lead); validation (equal); writing -original draft (supporting); writing -review and editing (equal).

ACK N OWLED G M ENTS
The authors thank the Farm Bill project. The authors also thank Prof.
Andrea Watzinger and DI. Katharina Schott for their invaluable guidance and assistance in the isotope laboratory, as well as Dr. Hanna Nadel and the USDA insectary for rearing the moth control groups.
Many thanks also go to Dr. Barbara Gepp, Johanna Wachter, and the OeAD for making the cooperation with the Reykjavik University possible. The main author also wishes to express her special thanks to Prof. Tommi Ekholm and Dr. Lukas Kohl for manuscript writing and data analysis advice.

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interest to disclose.

DATA AVA I L A B I L I T Y S TAT E M E N T
The moth and eggs metadata, measurement results and data analysis can be found under the DOI https://doi.org/10.5061/dryad.15dv4 1p0w.